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ABSTRACT 

We describe a new mechanism that leads to the destabilisation of non-axisymmetric 
waves in astrophysical discs with an imposed radial temperature gradient. This might 
apply, for example, to the outer parts of protoplanetary discs. We use linear density 
wave theory to show that non-axisymmetric perturbations generally do not conserve 
their angular momentum in the presence of a forced temperature gradient. This implies 
an exchange of angular momentum between linear perturbations and the background 
disc. In particular, when the disturbance is a low-frequency trailing wave and the 
disc temperature decreases outwards, this interaction is unstable and leads to the 
growth of the wave. We demonstrate this phenomenon through numerical hydrody¬ 
namic simulations of locally isothermal discs in 2D using the FARGO code and in 3D 
with the ZEUS-MP and PLUTO codes. We consider radially structured discs with a 
self-gravitating region which remains stable in the absence of a temperature gradient. 
However, when a temperature gradient is imposed we observe exponential growth of a 
one-armed spiral mode (azimuthal wavenumber m = 1) with co-rotation radius outside 
the bulk of the spiral arm, resulting in a nearly-stationary one-armed spiral pattern. 
The development of this one-armed spiral does not require the movement of the cen¬ 
tral star, as found in previous studies. Because destabilisation by a forced temperature 
gradient does not explicitly require disc self-gravity, we suggest this mechanism may 
also affect low-frequency one-armed oscillations in non-self-gravitating discs. 

Key words: accretion, accretion discs, hydrodynamics, instabilities, methods: nu¬ 
merical, protoplanetary discs 


1 INTRODUCTION 


An exciting development in the study of circumstellar discs 
is the direct observation of large-scale, non-axisymmetric 
structures within them. These include lopsided du st 


201a 

: iGasassus et al.l l2013l: llsella et al.l l2013l: Perez et al. 

2014 

: Follette et al.l 

2OI4I; van der Plas et al.l 2014 

\ and 

spiral arms (iHashimoto et al.l 1201 ll: iMuto et al.l 

2OI2I; 

Boccaletti et al.ll20ia 

: iGradv et al. 20131: Ghristiaens et al. 


20141: Avenhaus et al.ii2014). 


The attractive explanation for asymmetries in circum¬ 
stellar discs is disc-planet interaction. In particular, spiral 
structures naturally arise from the gravitational interaction 
between a planet and t he gaseous protoplane tary disc it is 
embedded in (see, e.g. iBaruteau et al.l l2013l . for a recent 
review). Thus, the presence of spiral arms in circumstel- 
lar discs could be signposts of planet formation (but see 
ljuhasz et al.l[2014l ). 
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Spiral arms are also characteristic of global gravi- 


(iGoldreich & Lvnden-Bellll965 

:lLaughlin & Rozvczkall99€ : 

Laughlin. Korchagin & Adams 

19981: 

Nelson et al. Il998: 

Lodato & Ricdl2005l: iForgan et al.ll2011 

). Large-scale spira’ 


arms can provide significant angul ar momentum transport 
neces s ary for mass accret i on |Lvnden-Bell fc Kalnais 
I972I: Papaloizou k, Savoniid Il99ll : iBalbus fc PapaloizOT 


I999I : iLodato fc Ricd 2004h . and spiral structures due to 


GI are potentially observable w ith the Atacama Large 
Milli m eter/sub-millim e ter A rray JCossins. Lodato fc Testa 
I2OI0I : iDipierro et al.l l2014l ) . GI can be expected in 
the earliest st ages of circumstellar disc formation 


iKratter et al.l 20101: Inutsuka. Machida fc Matsumotol 


I2OIOI: iTsukamoto. Machida fc Inutsuk^ l2013l). and ma; 


iga inutsuKal izuiai. ana may 
be possible in the outer parts of the disc (IRafikovI l2005l : 
iMatzner &: Levinll2005l : iKimura &: Tsuribell2012| ). 


Single-arm spirals, or eccentric modes, corresponding 
to perturbations with azimuthal wavenumber m = 1, have 
received interest in the context of protoplanetary discs 
because of their global nature (I Adams. Ruden fc Shul 
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198^ _ Heem sker k^ P apaloizou fc Savoniid 

Laughlin fc KorchaginI [l996l : ITi 


, _ .19921: 

_ __ __ _in d I 2 OO 1 I : IPapaloiz^ 

200d: Hopkinsl l201flll . In the ‘SLING’ mechanism proposed 
by IShu et al. lll990ll . an m = 1 gravitational instabil¬ 
ity arises from the motion of the central star induced 
by the one-armed perturbation, and requires a massive 
disc (the former may have observable consequences, 
iMichael fc Durisenll201CTl . 


In this work we identify a new mechanism that leads 
to the growth of one-armed spirals in astrophysical discs. 
We show that when the disc temperature is prescribed 
(called locally isothermal discs), the usual statement of 
the conservation of angular momentum for linear pertur¬ 
bations acquires a source term proportional to the tem¬ 
perature gradient. This permits angular momentum ex¬ 
change between linear perturbations and the background 
disc. This ‘background torque’ can destabilise low-frequency 
non-axisymmetric trailing waves when the disc temperature 
decreases outwards. 

We employ direct hydrodynamic simulations using three 
different grid-based codes to demonstrate how this ‘back¬ 
ground torque’ can lead to the growth of one-armed spirals 
in radially structured, self-gravitating discs. This is despite 
the fact that our disc models do not meet the requirements 
for the ‘SLING’ mechanism. Although our numerical simula¬ 
tions consider self-gravitating discs, this ‘background torque’ 
is generic for locally isothermal discs and its existence does 
not require disc self-gravity. Thus, the destabilisation effect 
we describe should also be applicable to non-self-gravitating 
discs. 

This paper is organised as follows. In Sj2] we describe 
the system of interest and list the governing equations. In 
j}51 we use linear theory to show how a fixed temperature 
profile can destabilise non-axisymmetric waves in discs. !j4] 
describes the numerical setup and hydrodynamic codes we 
use to demonstrate the growth of one-armed spirals due to 
an imposed temperature gradient. Our simulation results are 
presented in 35]for two-dimensional (2D) discs and in jjUlfor 
three-dimensional (3D) discs, and we further discuss them 
in m We summarise in !j8]with some speculations for future 
work. 


2 GOVERNING EQUATIONS 

We consider an inviscid fluid disc of mass Md orbiting a 
central star of mass M,. We will mainly examine 2D (or 
razor-thin) discs in favour of numerical resolution, but have 
also carried out some 3D disc simulations. Hence, for gen¬ 
erality we describe the system in 3D, using both cylindrical 
{R, 4>, z) and spherical polar coordinates (r, 9, (f>) centred on 
the star. The governing fluid equations in 3D are 

+ V ■ {pv) = 0, (1) 

^ + v ■ Vu ^ —-\/p — V<E>tot, (2) 

ot p 

V"<E>d = 47rGp, (3) 

where p is the mass density, v is the fluid velocity (the angu¬ 
lar velocity being Q = v^/R), p is the pressure and the total 
potential "Ltot = *!>* -I- 'Ld consists of that from the central 


star. 


'I>*(r) = — 


GM* 


(4) 


where G is the gravitational constant, and the disc potential 
$d. We impose a locally isothermal equation of state 

p = CsiR)p, (5) 


where the sound-speed Cs is given by 


Cs (R) — CsO 



-q/2 


( 6 ) 


where Cso = hRo^lk{Ro) and h is the disc aspect-ratio at 
the reference radius R = Ro, 0,k = \/GM^/R^ is the mid¬ 
plane Keplerian frequency, and —q is the imposed tempera¬ 
ture gradient since, for an ideal gas the temperature is pro¬ 
portional to Cj. For convenience we will refer to as the 
disc temperature. The vertical disc scale-height is defined 
by iL = Cs/Dfc. Thus, a strictly isothermal disc with g = 0 
has H oc R?^^, and q = 1 corresponds to a disc with constant 
H/R. 

The 2D disc equations are obtained by replacing p 
with the surface mass density S, p becomes the vertically- 
integrated pressure, and the 2D fluid velocity v is evaluated 
at the midplane, as are the forces in the momentum equa¬ 
tions. In the Poisson equation, p is replaced by 'E5{z), where 
S{z) is the delta function. 


3 LINEAR DENSITY WAVES 

We describe a key feature of locally isothermal discs that 
enables angular momentum exchange between small distur¬ 
bances and the background disc through an imposed radial 
temperature gradient. This conclusion results from the con¬ 
sideration of angular momentum conservation within the 
framework of linear perturbation theory. For simplicity, in 
this section we consider 2D discs. 

In a linear analysis, one assumes a steady axisymmetric 
background state, which is then perturbed such that 

E —>• E(R)-I-(5E„i(R) exp [i (—crt-I-mi)))], (7) 

and similarly for other variables, where cr = cu -|- iq is a 
complex frequency with w being the real frequency, 7 the 
growth rate, and m is an integer. We take m > 0 without 
loss of generality. 

The linearised mass and momentum equations are 

- ia^E^ {R^SvHrn ) - ^ , (8) 

- iaSvRm - 2mv^m = -g(R)^ (9) 

-iaSv^rn + ^SVRm = +S^rr?j , ( 10 ) 

where a = a — niQ and — R~^dR{R‘^V?) is the square 
of the epicyclic frequency. A locally isothermal equation of 
state has been assumed in Eq. [51 The linearised Poisson 
equation is 

V^i5-1>^ = 47rG(IE^5(z). (11) 

These linearised equations can be combined into a single 
integro-differential equation eigenvalue problem. We defer a 



























full numerical exploration of the linear problem to a future 
study. Here, we discuss some general properties of the linear 
perturbations. 


3.1 Global angular momentum conservation for 
linear perturbations 


It can be shown that linear perturbations with ^-dependence 
in the form exp (imc^) satisfies angular momentum conser¬ 
vation in the form 


djw 


dt 

(e-K- 


+ V • = Tbg, 


Naravan. Gold reich &: GoodmanI 


mS . 


di 


ilin =- — Im ( — -I- nfc • ^ X -I- imfll^l" 


( 12 ) 

19871: 


IRvu fc GoodmanI Il992l : Ibin. Panaloizou fc Kiev! Il99 
where 


(13) 


is the angular momentum density of the linear disturbance 
(which may be positive or negative), ^ is the Lagrangian dis¬ 
placement and * denotes complex conjugation, and F is the 
vertically-integrated angular momentum flux consisting of a 
Reyn olds stress and a gravitational stress llLin &: Panaloizoul 
201l|) . Explicit expressions for ^ can be found in, e.g. 
Papaloizou fc Pringld lll985l l. 

In Eq. [H the background torque density Tbg is 

= . (14) 


and arises because we have adopted a locally isothermal 
equation of state in the perturbed disc. We outline the 
derivation of Tbg in Appendix m 

In a barotropic fluid, such as a strictly isothermal disc, 
Tbg vanishes and the total angular momentum associated 
with the perturbation is conserved, provided that there 
is no net angular mom entum flux. However, as noted in 
Ibin fc Papaloizo'ul ll201lll . if there is an imposed temperature 
gradient, as in the disc models we consider, then Tbg 0 in 
general, which corresponds to a local torque exerted by the 
background disc on the perturbation. 

The important consequence of the background torque is 
the possibility of instability if TuGiun > 0. That is, if jun is 
positive (negative) and Tbg is also positive (negative), then 
the local angular momentum density of the linear distur¬ 
bance will further increase (decrease) with time. This implies 
the amplitude of the disturbance may grow by exchanging 
angular momentum with the background disc. 

We demonstrate instability for low-frequency modes 
(|tj| <C mfl) by explicitly showing its angular momentum 
density jn-n < 0, and the background torque Tbg < 0 for ap¬ 
propriate perturbations and radial temperature gradients. 


3.2 Angular momentum density of 

non-axisymmetric low-freqnency modes 

From Eq. and assuming a time-dependence of the form 
exp (—iut) with Reu = lj, the angular momentum density 
associated with linear waves is 

jun = ^[(c^- mn) 1 ^ 1 = + iff (e^e; - Cr^^)] ■ (15) 
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For a low-frequency mode, |a;| mfl. Then neglecting the 
term we And 

iiin ~ -f i (Cfle; - Cr^ci,)] 

= [("i-1)1^1"+ IC« + iec^l"] ■ (16) 

Thus, non-axisymmetric (m ^1) low-frequency modes 
generally have negative angular momentum. If the mode 
loses (positive) angular momentum to the background, then 
we can expect instability. We show below how this is possible 
through a forced temperature gradient via the background 
torque. It is simplest to calculate Tbg in the local approxi¬ 
mation, which we review first. 


3.3 Local results 


In the local approximation, perturbations are assumed to 
vary rapidly relative to any background gradients. The dis¬ 
persion relation for tightly-wound density waves of the form 
exp [i(—erf -|- mi^ -f kR)] in a razor-thin disc is 

{a-mnf = k'^+ Oc^,-2TvGT:\k\, (17) 

wher e fc is a real wavenumber such that \kR\ 3> 1 (IShul 
Note that in the strictly local approximation, where 
all global effects are neglected, only axisymmetric perturba¬ 
tions (m = 0) can be unstable. 

Given the real frequency uj or pattern speed Q.p = cu/m 
of a non-axisymmetric neutral mode , Eq. [T7] can be solved 
for |fe|. 


\k\ = fc. [l ± Vi - Q"(l - i^")] , 

where 
. ttGE 


is a characteristic wavenumber, 


Q 


CsR 

ttGE 


(18) 

(19) 

( 20 ) 


is the usual Toomre parameter, and 


(w — m^l) 

K 


( 21 ) 


is a dimensionless frequency. In Eq. 1181 the upper (lower) 
sign correspond to short (long) waves, and A: > 0 (A: < 0) 
correspond to trailing (leading) waves. 

At the co-rotation radius Rc the pattern speed matches 
the fluid rotation. 


Q{Rc) = flp. 


( 22 ) 


Lindblad resonances Rl occurs where 

v\Rl) = 1. (23) 


Finally, Q-barriers occur at radii Rqs where 

Q^iRQi) [1 - ^^RQi)] = 1. (24) 

According to Eq. 1181 purely wave-like solutions with real k 
are only possible where Q^{1 — ^ 1. 

A detailed di scuss i on of the properties of local density 
waves is given in IShul lll99lll . An important result, which 
holds for waves of all frequencies, is that waves interior to co¬ 
rotation [R < Rc) have negative angular momentum, while 
waves outside co-rotation [R > Rc) have positive angular 
momentum. 


































4 


Lin 


3.4 Unstable interaction between low-frequency 
modes and the background disc due to an 
imposed temperature gradient 

Here we show that the torque density acting on a local mode 
due to the background temperature gradient can be neg¬ 
ative, which would enforce low-frequency modes, because 
they have negative angular momentum. 

The Eulerian surface density perturbation is given by 

= = (25) 

We invoke local theory by setting d/dR —>■ ifc where k is real, 
and assume \kR\ ';$> m so that the second term on the right 
hand side of Eq. [25] can be neglected. Then 

STim — —ifcE^fl. (26) 


Inserting this into the expression for the background torque, 
Eq. 1141 we find 


„ m dc?. , ,o 


(27) 


This torque density is negative for trailing waves {k > 0) 
in discs with a fixed temperature profile decreasing outwards 
{dcl/dR < 0). Note that this conclusion does not rely on the 
low-frequency approximation. 

However, if the linear disturbance under consideration 
is a low-frequency mode, then it has negative angular mo¬ 
mentum. If it is trailing and dc^/dR < 0, as is typical in 
astrophysical discs, then Tbg < 0 and the background disc 
applies a negative torque on the disturbance, which further 
decreases its angular momentum. This suggests the mode 
amplitude will grow. 

Using jiin and Tbg, we can estimate the growth rate 7 
of linear perturbations due to the background torque as 


27 ~ 


Tbg 

jlin 


(28) 


where the factor of two accounts for the fact that the angular 
momentum density is quadratic in the linear perturbations. 
Inserting the above expressions for jun and Tbg for gives 


dcj k _ IC-Rp _ 

di?H[(m-l)|^|2 + |^« + ie,^|2] 


dc^ k 
dR mil ’ 


(29) 


where the second equality uses ~ 2i^jijm for low- 
frequency modes in the local approximation, as shown 
in Appendix E Then for the temperature profiles = 
Cso{R/Ro)~'^ as adopted in our disc models, 


where we used Ca ~ Hil ~ hRil. Ea. l30l suggests that pertur¬ 
bations with small radial length-scales {kH > m) are most 
favourable for destabilisation. Taking the local approxima¬ 
tion is then appropriate. 

Note that the derivation of Eq. [27] and Eq. [30] do not 
require the disc to be self-gravitating. Thus, destabilisation 
by the background torque is not directly associated with disc 
self-gravity. However, in order to evaluate Eq. [23 or Eq. [30] 
in terms of disc parameters (as done in 45.311 . we need to 
insert a value of k, which may depend on disc self-gravity 
(e.g. from Eg. [TSll. 


4 NUMERICAL SIMULATIONS 

We demonstrate the destabilising effect of a fixed tempera¬ 
ture gradient using numerical simulations. The above discus¬ 
sion is generic for low-frequency non-axisymmetric modes, 
but for simulations we will consider specific examples. 

There are two parts to the destabilising mechanism: 
the disc should support low-frequency modes, which is then 
destabilised by an imposed temperature gradient. The lat¬ 
ter is straight forward to implement by adopting a locally 
isothermal equation of state as described in ^ For the for¬ 
mer, we consider discs with a radially-structured Toomre Q 
profile. We use local theory to show that such discs can trap 
low-frequency one-armed (m = 1) modes. This is convenient 
because Eq. [30] indicates that modes with small m are more 
favourable for destabilisation. 


4.1 Disc model and initial conditions 


For the initial disc profile we adopt a modified power-law 
disc with surface density given by 


E(R) = E„f X B{R-, Ri,R2, e, SR), (31) 

where s is the power-law index describing the smooth disc 
and Eref is a surface density scale chosen by specifying Qout, 
the Keplerian Toomre parameter at R = R 2 , 


Qout - 


Csilk 


ttGE 


R—R 2 


(32) 


The bump function B{R) represents a surface density boost 
between R £ [-Ri,-R 2 ] by a factor > 1, and SR is the 
transition width between the bump and the smooth disc. 
We choose 


B(R) = /i(R) x/2(R), 


HR) = 5 (1 - ^) 


1 -I- tanh 
1 — tanh 


(^) 

(^) 


+ 

+ 


(33) 

(34) 

(35) 


where Ai ,2 = SR x i 7 (-Ri, 2 ). 

The 3D disc structure is obtained by assuming vertical 
hydrostatic balance 


I dp d4»« 

p dz dz dz 


(36) 


which gives the mass density as 




(37) 


where Z{R,z) describes vertical stratification. In practice, 
we numerically solve for Z{R,z) by neglecting the radial 
self-gravity force compared to vertical self-gravity, which re¬ 
duces the equations for vertical hydrostatic equilibrium to 
or dina r y diff erential equations. This procedure is described 
inlLhil J2012lb 


Our fiducial parameter values are: s = 2, Ri = Rq, 
R 2 = 2Ro, e = 0.1, SR — 5, h — 0.05 and Qout = 2. An 
example of the initial surface density and the Toomre Q 
parameter is shown in Fig. [T] Since Q > 1, the disc is sta¬ 
ble to local axisymmetric perturbations llToomrdll964h . The 
transition between self-gravitating and non-self-gravitating 
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Figure 1. Fiducial profiles of the surface density (top) and 
Toomre parameter (bottom) used in this work. 


portions of the disc occur smoothly across ~ lOH. Initially 
there is no vertical or radial velocity {vn = Vr = vg = 0). 
The azimuthal velocity is initialized to satisfy centrifugal 
balance with pressure and gravity, 

r p or or 
and similarly in 2D. 


4.2 Codes 

We use three independent grid-based codes to simulate the 
above system. We adopt computational units such that G — 
Mt = Rq = 1. Time is measured in the Keplerian orbital 
period at the reference radius, Po = 27r/Dj;(i?o)- 


The 2D Poisson equation is solved in integral form, 


^d,z=o{R, 4’) 



G'E{R',cj)')R'dR'd(t>' 

R? + R'^ — 2RR' cos {(j) — (f>') + ’ 


(39) 


using Fast Fourier Transform (FFT), where tg is a softening 
length to prevent a numerical singularity. Th e FFT approach 
requires Eg oc R (iBaruteau fc Massedl2008h . In FARGO, Eg 
is set to a fraction of hR. 


4.2.S ZEUS-MP 

ZEUS-MP is a general-purpose finite difference code 
llHaves et al.ll200^ ). We use the code in 3D spherical ge¬ 
ometry, covering r £ [rmin, rmax], 0 £ [6»niin, vr/2], <!> £ [0, 27r]. 
The vertical domain is chosen to cover hh scale-heights at 
R = Ro, i.e. tan (7r/2 — 0min)/h = uh - The grid is loga¬ 
rithmically spaced in radius and uniformly spaced in the 
angular coordinates. We assume symmetry across the mid¬ 
plane, and apply reflective boundary conditions at radial 
boundaries and the upper disc boundary. 

ZEUS-MP solves the 3D Poisson equation using a conju¬ 
gate gradient method. To supply boundary conditions to the 
linear solver, we expand the bo undary pote ntial in spherical 
harmonics Yim as described in IBos i (1 19801) . The expansion 
is truncat ed at {l,m) = (lmax,mmax). This code was used in 
iLinl (|2012l ) for self-gravitating disc-planet simulations. 


4.2.3 PLUTO 

PLU TO is a general-purpose Godunov code dMignone et al.l 
l2007h . The grid setup is the same as that adopted in ZEUS- 
MP above . We configure the code similarly to that used in 
[Li^(l2014l) : piece-wise linear reconstruction, a Roe solver and 
second order Runge-Kutta time integration. We also enable 
the FARGO algorithm for azimuthal transport. 

We solve the 3D Poisson equation thr oughout the do¬ 
main using spherical harmonic expansion IIBos , as 

used for the boundary pot ential in ZEUS - MP. T his ver¬ 
sion of PLUTO was used in iLin fc Gloutierj (l2014h for self- 
gravitating disc-planet simulations, producing similar re¬ 
sults to that of ZEUS-MP and PARGO. 


4.2.1 FARGO 

Our primary cod e _ is FARGO with self-gravity 

(IBaruteau fc Mass^ l2008l ) . This is a popular, simple 
finite-difference code for 2D discs. ‘FARGO’ refers to its 
azimuthal transport algorithm, which removes the mean 
azimuthal velocity of the disc, thereby permit larger time- 
steps than that would otherwise be allowed by the usual 
Courant cond ition based on the full azimuthal velocity 
(^ asse ^ l2000al l^. 

The 2D disc occupies R £ Rmax], <?!> £ [0, 27r] and 

is divided into (Nu, A^</>) grids, logarithmically spaced in ra¬ 
dius and uniformly spaced in azimuth. At radial boundaries 
we set the hydrodynamic variables to their initial values. 


4.3 Diagnostics 

4 . 3.1 Evolution of non-axisymmetric modes 

The disc evolution is quantified using mode amplitudes and 
angular momenta as follows. We list the 2D definitions with 
obvious 3D generalisations. A hydrodynamic variable / (e.g. 
E) is written as 


00 

f{R,(j),t)= ^ /m {R, t) exp im<j) 

m= — oo 


— /o + 2 Re 


fm exp (im<?!>) , 


(40) 


_m = l 

where the fm may be obtained from Fourier transform in (f. 
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The normalised surface density with azimuthal 
wavenumber m is 

AE^ = Re [Em exp (im0)] (41) 

Eoo 

where Eoo = ^o{t = 0). The time evolution of the m**' 
mode can be characterized by assuming Em oc exp {—iat) as 
in linear theory. The total non-axisymmetric surface density 
is 

AE = (42) 



t/P„ 


4-3.2 Angular momentum decomposition 


The total disc angular momentum is 
J = / / 'ERvcf)RdRd(j) 

/ -Rmax 

RT,QV^oRdR 

^min 

oo /. R CO 

/ -‘‘■max 

^ 27r / 2RRe [Y^mV^rn] RdR = ^ (43) 

m = l ^min m = 0 




We will refer to Jm as the m**' component of the total an¬ 
gular momentum, and use it to monitor numerical angular 
momentum conservation in the simulations. It is important 
to distinguish this empirical definition from the angular mo¬ 
mentum of linear perturbations given in ^ which is defined 
through a conservation law. 


4-3.3 Three-dimensionality 


In 3D simulations we measure the importance of vertical 
motion with 0, where 


_ {vl 

~ K) + {viy 


(44) 


and (•) denotes the density-weighted average, e.g., 

{vL = 




r ^2 r 


J^PdV 


( 45 ) 


and similarly for the horizontal velocities. Thus 0 is the 
ratio of the average kinetic energy associated with vertical 
motion to that in horizontal motion. The radial range of 
integration is taken over r £ [Ri,R 2 \ since this is where we 
find the perturbations to be confined. 


5 RESULTS 

We first present results from FARGO simulations. The 2D 
disc spans [Rmin, Rmax] = [0.4,10]Ro- This gives a total disc 
mass Md = 0.086M*. The mass within R £ [i?min, Ri] is 
0.017M*, that within R £ [Ri,R 2 ] is 0.049M*, and that 
within R £ [R 2 ,Rmax] is 0.021M*. We use a resolution of 
Nr X = 1024 X 2048, or about 16 grids per R^nd adopt 
eg — \Q~*hR for the self-gravity softening lengtly. 


^ In 2D self-gravity, eg also approximates for the vertical disc 
thickness, so a more appropriate value would be Cg ~ R 


Figure 2. Evolution of non-axisymmetric surface density maxima 
in the FARGO simulation initialised with perturbations with m € 
[ 1 , 10 ). 


In these simulations the disc is subject to initial pertur¬ 
bations in cylindrical radial velocity, 

M 

COS mi)), ( 46 ) 

m=l 

where the amplitude 5 £ [—10“®, 10“®] is set randomly but 
independent of 0, R = (Ri -I- R2)/2 and AR = (R 2 — Ri)/2. 


VR ^ VR-\- Ca — exp 


1 fR-R 
'2 V AR 


5.1 Reference run 


To obtain a picture of the overall disc evolution, we de¬ 
scribe a fiducial run initialised with M = 10 in Eq. 1461 Fig. 
[ 2 ] plots evolution of the maximum non-axisymmetric sur¬ 
face density amplitudes in R £ [Ri,R 2 ] for m £ [1,10]. 
Snapshots from the simulation are shown in Fig. [3] At 
early times t < lOORo the disc is dominated by low- 
amplitude high-m perturbations. The m ^ 4 modes growth 
initially and saturate (or decays) after t = 40Ro. No¬ 
tice the low m ^ 2 modes decay initially, but grows be¬ 
tween t £ [20,40]Ro, possibl y due to non-linear interac¬ 
tion of the high-m modes IILaughlin fc Korchagin! Il996l : 
iLaughlin, Korchagin fc Adams! 1 1997 1. However, the m = 1 
mode begins to grow again after t = 70Ro, and eventually 
dominates the annulus. 

Fig. [4] shows the evolution of disc angular momentum 
components. Only the m = 0, 1 components are plotted 
since they are dominant. The m = 1 structure has an as¬ 
sociated negative angular momentum, which indicates it is 
a low-frequency mode. Its growth is compensated by an in¬ 
crease in the axisymmetric component of angular momen¬ 
tum, such that A Jo -I- A Ji ~ 0. Note that FARGO does not 
conserve angular momentum exactly. However, we find the 
total angular momentum varies by |AJ/J| = 0(10“®), and 
is much smaller than the change in the angular momenta 
components, |AJo,i/J| > 0(10“®). Fig. |4]then suggest that 
angular momentum is transferred from the one-armed spiral 
to the background disc. 


dMuller. Kiev fc Merul[2012h . However, because eg oc R is needed 
in FARGO, the Poisson kernel (Eq. I39II is no longer symmetric 
in {R, R'). We choose a small eg in favour of angular momentum 
conservation, keeping in mind that the strength of self-gravity 
will be over-estimated. 
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Figure 3. Visualisation of the FARGO 2D simulation in Fig. [2] The total non-axisymmetric surface density AS is shown. 



Figure 4. Evolution of angular momentum components in the 
FARGO simulation in Fig. [2H3 The perturbation relative to t = 
0 in 2D is shown in units of the initial total angular momentum 
.^ref * 

5.2 Dependence on the imposed temperature 
profile 

We show that the growth of the m = 1 spiral is associ¬ 
ated with the imposed temperature gradient by performing 
a series of simulations with q £ [0,1]. However, to main¬ 
tain similar Toomre Q profiles, we adjust the surface density 
power-law index such that s = (3 -I- q)/2. For clarity these 
simulations are initialised with m = 1 perturbations only. 

Fig. [5] compares the m = 1 spiral amplitudes as a func¬ 
tion of q. We indeed observe slower growth with decreas¬ 
ing q. Although the figure indicates growth for the strictly 
isothermal disc (g = 0 ), we did not observe a coherent one- 
armed spiral upon inspection of the m = 1 surface density 
held. The growth in this case may be associated with high-m 
modes, which dominated the simulation. 

We plot growth rates of the m = 1 mode as a function 
of q in Fig. |6l The correlation can be htted with a linear 
relation 

7 ~ [0.015q - 7.9 X lO”"*] flkiRo)- 


m=1 



Figure 5. Evolution of the m = 1 spiral amplitude as a func¬ 
tion of the negative of the imposed temperature gradient q. The 
maximum value of the m = 1 surface density in A S [iJi, R 2 ] is 
shown. 


As the background torque is proportional to q (Eq. I30II , this 
indicates that the temperature gradient is responsible for 
the development of the one-armed spirals observed in our 
simulations. 

We also performed a series of simulations with variable 
aspect-ratio h £ [0.03, 0.07] but hxed q = 1. This affects 
the magnitude of the temperature gradient since Cs oc h. 
However, with other parameters equal to that in the fiducial 
simulation, varying h also changes the disc mass. For h £ 
[0.03, 0.07] the total disc mass ranges from Md = 0.052M* 
to Md = 0.12M* and the mass within R £ [J?i,i? 2 ] ranges 
from 0.033M* to 0.062M*. 

Fig. (3 shows the growth rates of the m = 1 spiral in 
R € [i?i, R 2 \ as a function of h. Growth rates increases with 
h, roughly as 

7 ~ [O.lO/i -b 8.3 X 10"®] nk{Ro)- 

However, a linear ht is less good than for variable q cases 
above. This may be due to the change in the total disc mass 
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Figure 6. Growth rates of the m = 1 spiral mode as a function 
of the imposed sound-speed gradient q (asterisks). A linear fit is 
also plotted (dotted line). 
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Figure 7. Growth rates of the m = 1 spiral mode as a function 
of the disc aspect-ratio h (asterisks). A linear fit is also plotted 
(dotted line). 
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when h changes. We find no qualitative difference between 
the spiral pattern that emerges. 



Figure 8. Cartesian visualisation of the m = 1 surface density 
structure in the FARGO simulation initialised with only m = 1 
perturbations. 



5.3 Properties of the m = 1 spiral and its growth 

Here we analyse the 9 = 1 case in Fig. [S]in more detail. 

Fig. [ 8 ] shows a snapshot of the m = 1 surface density of 
this run. By measuring the m — 1 surface density amplitude 
and its pattern speed, we obtain a co-rotation radius and 
growth rate 


Rc — 4.4iZo, 

7 ~ 0.0140fc(iZo) = O.130p. 


This one-armed spiral can be considered as low frequency 
because its pattern speed Op ~ 0.10fc(i?o) 0.30 in 

R G [i?i, R 2 ] (where it has the largest amplitude). Thus, the 
spiral pattern appears nearly stationary. The growth rate 
7 is also slow relative to the local rotation, although the 
characteristic growth time 7 “^ ~ lOPo is not very long. 

Next, we write Si = |Ei| exp (ifcP), where k is real, and 
assume the amplitude |Ei| varies slowly compared to the 
complex phase. This is the main assumption in local theory. 
We calculate k numerically and plot its normalised value in 
Fig. [9l We find 


kR 


1 

ci hQ ’ 


Figure 9. Normalised radial wavenumber of the m = 1 spiral in 

Fig. El 


where we used Q ~ CsO./nG'E and R£l/cs ~ h~^. Since 
Q = 0(1) and h <C 1 imply \kR\ ';$> 1, we can apply results 
from local theory ( 1)3.311 . Note also that k ~ k^. Fig.| 8 ]shows 
the m — 1 spiral is trailing, consistent with k > 0. 

Using the estimated value of Rc, we plot in Fig. llOl the 
quantity — 1 -\- Q~^, which is required to be positive in 
local theory for purely wave-like solutions to the dispersion 
relation (Eq. I17II when the mode frequency is given. Fig. HOI 
shows two Q-barriers located in the inner disc, at Rqi, = Ro 
and Pq(, = 1.6Po; the bounded region is indeed where the 
m — 1 spiral develops. This shows that the one-armed spiral 
is trapped. Note in this region, — ~ 0.1 <C 1, which 

is necessary for consistency with the measured wavenumber 
k and Eg. 1181 There is one outer Lindblad resonance at P_l ~ 
7.2Ro. Thus, acoustic waves may be laun ched in R > 7.2Ro 
by th e spiral disturbance in the inner disc IILin fc Papaloizoul 

We can estimate the expected growth rate of the m = 1 
mode due to the temperature gradient. Setting k = kc and 
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Figure 10. Dimensionless mode frequency u for the m = 1 spiral 
in Fig. [8] For a given real mode frequency, the dispersion relation 
for local density waves, Eq. [m permits purely wave-like solutions 
in regions where — 1 Q~^ > 0. 


m=1 



Figure 11. Rate of change of the m = 1 wave angular momentum 
as defined by Eq. 1481 (solid) compared to the torque exerted on 
the wave associated with the background temperature gradient 
(dotted). 


m = 1 into Eq. [30] gives 


7 



(47) 


Inserting q = 1, h = 0.05 and Q ~ 1.5 gives 7 ~ 0.01717, 
consistent with numerical results. 


5.3.1 Angular momentum exchange with the background 
disc 

We explicitly show that the growth of the m = 1 spiral is 
due to the forced temperature gradient via the background 
torque described in ^Sl We integrate the statement for an¬ 
gular momentum conservation for linear perturbations, Eq. 
n assuming boundary fluxes are negligible, to obtain 

J rRmt^x rRmax 

— / jii^2nRdR^ / TnG^TvRdR, (48) 

7 ^min "I 

4lin 

where we recall Tbg is the torque density associated with 
the imposed sound-speed profile (Eq. I14II . We compute both 
sides of Eq. [48] using simulation data, and compare them in 
Fig. El There is a good match between the two torques, es¬ 
pecially at early times t < llOPo- The average discrepancy is 
~ 5%. The match is less good later on, when the spiral am¬ 
plitude is no longer small (maxAEi ~ 0.2 a.t t = llOPo and 
maxAEi ~ 0.4 hy t = 120Po) and linear theory becomes 
less applicable. Fig. El confirms that the m = 1 spiral wave 
experiences a negative torque that further reduces its (neg¬ 
ative) angular momentum, leading to its amplitude growth. 
This is consistent with angular momentum component mea¬ 
surements (Fig. 3)). 


6 THREE-DIMENSIONAL SIMULATIONS 

In this section we review 3D simulations carried out using 
ZEUS-MP and PLUTO. The main purpose is to verify the 
above results with different numerical codes, and validate 
the 2D approximation. 

The 3D disc has radial size [rmin, ^max] = [0.4, 10]Po 
and vertical extent uh = 2 scale-heights at i? = Pq. The 
resolution is Nr x Ne x = 256 x 32 x 512, or about 


4 cells per H. Because of the reduced resolution compared 
to 2D, we use a smooth perturbation by setting 5 = 10~® 
and M = 1 in Eq. 1461 This corresponds to a single m = 1 
disturbance in P € [Pi, P 2 ]. 

The 3D discs are initialised in approximate equilib¬ 
rium only, so we first evolve the disc without perturba¬ 
tions using (imax, m-max) = (32,0) Up to t = IOPq, dur¬ 
ing which meridional velocities are damped out. We then 
restart the simulation with the above perturbation and 
(^max, nimax) = (32, 32), and damp meridional velocities near 
the radial boundaries. 

Fig. El plots the evolution of the m = 1 spiral ampli¬ 
tudes measured in the ZEUS-MP and PLUTO runs. We also 
ran simulations with a strictly isothermal equation of state 
(q = 0), which display no growth compared to that with a 
temperature gradient. This confirms the temperature gradi¬ 
ent effect is the same in 3D. 

In the ZEUS-MP run, we observed high-m disturbances 
developed near the inner boundary initially, which is likely 
responsible for the growth seen at t < 50Po. This is a nu¬ 
merical artifact and effectively seeds the simulation with a 
larger perturbation. Results from ZEUS-MP are therefore 
off-set from PLUTO by ~ 50Po. However, once the coherent 
m = 1 spiral begins to grow (t > lOOPo), we measure similar 
growth rates in both codes: 

7 ~ 0.0073Dfc (Po) PLUTO, 

7 ~ 0.0085Dfc (Po) ZEUS-MP. 

Both are somewhat smaller than the 2D simulations. This 
is possibly because of the lower resolutions adopted in 3D 
and /or because the effective Too mre parameter is larger in 
3D (iMamatsashvili fc RicdBoid ') which, from Eq. 1471 is sta¬ 
bilising. 

Visualisations of the 3D simulations are shown in Fig. 
for the disc midplane and near the upper disc boundary. 
The snapshots are chosen when the one-armed spirals in the 
two codes have reached comparable amplitudes. Both codes 
show similar one-armed patterns at either height, and the 
midplane snapshot is similar to the 2D simulation (Fig. [3]). 
The largest spiral amplitude is found in the self-gravitating 
region P G [Pi,P 2 ], independent of height. However, notice 
the spiral pattern extends into the non-self-gravitating outer 
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m=1 



Figure 12. Evolution of the maximum m = 1 density component 
in r E [Ri,R 2 ] in the 3D simulations. Results from discs with 
a temperature gradient (^q = 1) and a strictly isothermal disc 
= 0) are shown. 
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Figure 13. Three-dimensional simulations using the ZEUS-MP 
(top) and PLUTO (bottom) codes. The m = 1 density component 
Apr at the midplane (left) and approximately two scale-heights 
above the midplane (right) is shown. Here ip = 7r/2 — 0 is the 
angular height from the midplane. 


disc {R > R 2 ) at 2 ~ 2H, i.e. the disturbance becomes more 
global away from the midplane. 


6.1 Vertical structure 

Fig. [ 14 ] shows the vertical structure of the one-armed spi¬ 
ral in the PLUTO run. The spiral is vertically conhned to 
z ^ H at R Ro (the self-gravitating region). Thus, a 2D 
disc model, representing dynamics near the disc midplane, 
is sufficient capture the instability. However, for R > 2Rq 
the spiral amplitude increases away from the midplane. It 
remains small in our disc model (|Api| < 0.1), but could be¬ 
come signihcant with a larger vertical domain. This means 



Figure 14. The m = 1 density component in the merid¬ 
ional plane in the PLUTO simulation. The slices are taken at 
the azimuth of max[Api(r, 7r/2, (ji)]. Arrows represent the vector 
{vr/Ra,—vg/rhsir? 9). The top (bottom) panel corresponds to 
the inner (outer) portions of the disc. 


that 3D simulations are necessary to study the effect of the 
one-armed spiraf on the exterior disc. 

Afthough Fig. [14] appears to display signihcant verti¬ 
cal motion, we measure the three-dimensionality parameter 
0 < 10“^ (Ea. l44ll . so vertical motions are insignihcant com¬ 
pared to horizontal motions. This supports a 2D approxima¬ 
tion. On the other hand, we Hnd max\vz/cs\ ~ 0.2 which, 
although sub-sonic, is not very small. 


6.2 Angular momentum conservation 

Fig. [ 15 ] shows the angular momentum evolution in the 3D 
runs during the linear growth of the one-armed spiral. Be¬ 
cause the ZEUS-MP simulation is off-set from PLUTO, the 
time intervaf for the pfot was chosen such that the change in 
the angular momentum components are comparable in the 
two codes. 

ZEUS-MP does not conserve angular momentum very 
well, but the variation in total angular momentum | A J/J| < 
0(10“®) is small compared to the individual components 
|AJo,i/J| ~ 10“"^. PLUTO reaches similar values of |AJo,i|, 
but achieves better conservation, with |AJ/J| = 0(10“®). 
These plots are again similar to the 2D simulations, i.e. an¬ 
gular momentum lost by Ji is gained by Jo- This confirms 
that the interaction between Ji and Jo operates in 3D and 
2D similarly. 


7 DISCUSSION 

We discuss below several issues related to self-gravitating 
discs in the context of our numerical simulations. However, 
it is important to keep in mind that the growth of the one- 
armed spiral in our simulations is not a gravitational insta¬ 
bility in the sense that destabilisation is through the back- 








































One-armed spirals 11 



150 160 170 180 190 200 

t/Po 


Figure 15. Evolution of angular momentum components in the 
3D simulations. The perturbation relative to t = IOPqi during the 
growth of the one-armed spiral, is shown in units of the initial 
total angular momentum Jref- 


ground torque associated with a forced temperature gradi¬ 
ent, and not by self-gravitational torquefl 


7.1 Motion of the central star 


In our models we have purposefully neglected the in¬ 
direct potential associated with a non-inertial reference 
frame to avoid complications arising from the motion 
of the central star. Although it has been established 
that such m otion can destabilise an m = 1 disturbance 
in the disc (lAda rns. Ru den fc Shul 1 19891 : IShu et ahl Il99d : 
iMichael fc PurisCT 2 OI 0 II . the disc masses in our models 


{Md < O.IM,) are not expected to be sufficiently massive 
for this effect to be significant. Indeed, simulations including 
the indirect potential, carried out in the early stages of this 
project, produced similar results. 


7.2 Role of Llndblad and co-rotation torques 

One effect of self-gravity is to allow the one-armed spiral, 
confined to R ~ i?o in our models, to act as an external po¬ 
tential for the exterior disc i n R > Rp. This is analogou s 
to disc-satellite interaction (iGoldreich fc Tremaine! Il979ll . 
where the embedded satellite exerts a torque on the disc 
at Lindblad and co-rotation resonances. 

In Appendix 0 we estimate the magnitude of this ef- 
fect using basic result s from disc-planet theory (see, e.g. 
IPapaloizou et aLll2007l . and references therein). There, we 


^ In fact, additional simulations with Qout = 4 (giving Q > 2.5 
throughout the disc) still develops the one-armed spiral, but with 
a smaller growth rate. 


find that the angular momentum exchange between the one- 
armed spiral and the exterior disc is insignificant compared 
to the background torque. 

We confirmed this with additional FARGO simulations 
that exclude the co-rotation and outer Lindblad resonances 
(OLR) by reducing radial domain size to Rmax = 3Ro, which 
still developed the one-armed spiral. 


7.3 Applicability to protoplanetary discs 

7.3.1 Thermodynamic requirements 

A locally isothermal equation of state represents the ideal 
limit of infinitely short cooling and heating timescales, so 
the disc temperature instantly returns to its initial value 
when perturbed. The background torque is generally non¬ 
zero if the resulting temperature profile has a non-zero radial 
gradient. 

A short cooling timescale tc can occur in the outer 


Rice & Armitaee 

200d: iGossins. Lodato & Clarkd 

201 c; 

Tsukamoto et al. 

120151). However, if a disc with O - 

; 1 is 

cooled (towards 

ero temperature) on a timescale tc S. Hi. 

it will fragment following gravitational instabilitv (iGammie 

200 ll:lRice. Lodato & Armitagell2005l:IPaardekooDeill2012l). 


Fragmentation can be avoided if the disc is heated to 
maintain Q > Qc, th e threshold for frag mentation (Qc — 1-4 
for isothermal discs. iMaver et al]l2004l ). This may be pos¬ 
sible in the out er parts of protoplanetary discs due to stel¬ 


lar irradiation ( Rafikovll2009l : iKratter fc Murrav-Glavll201ll : 


IZhu et ani2012 ). Sufficiently strong external irradiation is 


expected to suppress the l inear gravitational instabilities al¬ 
together (iRice et al.ll201ll ). 

The background torque may thus exist in the outer 
parts of protoplanetary discs that are irradiated, such that 
the disc temp erature is set externally with a n on-zero radial 
gradient ('e.g. lStamatellos fc Whitworthll200^ . Of course, if 
external irr adiation sets a strictly isothermal outer disc (e.g. 
iBolevlIi^OOgl ). then the background torque vanishes. 


7.3.8 Radial disc structure 

In our simulations the m = 1 spiral is confined between 
two Q-barriers, where real solutions to the local dispersion 
relation is possible fEg. 1181) . The existence of such a cavity 
results from the adopted initial surface density bump (Eq. 
1331) . Thus, in our disc models the main role of disc structure 
and self-gravity is to allow a local m = 1 mode to be set up, 
which is then destabilised by the background torque. 

In order to confine an m = 1 mode between two Q- 
barriers, we should have — = 1 at two radii. Assum¬ 

ing Keplerian rotation and a slow pattern speed flp Q, 
this amounts to 

(^) =2Q\Rq,). (49) 

Then two Q-barriers may exist when the Q^ profile rises 
more rapidly (decays more slowly) than R~^^^ for decreasing 
(increasing) R. Note that Eq. [dadoes not necessarily require 
strong self-gravity if Rc is large. 

A surface density bump can develop in ‘dead zones’ 
of protoplanetary discs, where there is reduced mass 

































































12 Lin 


accretion because the magneto-rotational inst ability is 
ineffective for angu l ar momentum transpor t dGammid 
Il99fil : [Turner fc Sanol l2008l : iLandrv et alJ l201ij i. The dead 


zon e becomes self-gravitating with sufficient mass bui lt- 


UD ([Armitage. Livio & Pringle! 2001): Martin et al.l 2012a tJ: 

Zhu et al. 

2009[. [201fl[: [Zhu. Hartmann fc Gammid [20101: 

Bae et al. 

[2013f). 


However, conditions in a dead zone may not be suit¬ 
able for sustaining a background torque because it may not 
cool/heat fast enough to maint ain a fixed temperature pro¬ 
file. Recently, iBae et ^ (l2014l ~l presented numerical mod¬ 
els of dead zones including a range of heating and cool¬ 
ing processes, which show that dead zones developed large- 
scale (genuine) gravitational instabilities with multiple spi¬ 
ral arms. Although this does not prove absence of the back¬ 
ground torque, it is probably insignificant compared to grav¬ 
itational torques. 

Another possibility is a gap opened by an embedded 
planet. In that case Q rises rapidly towards the gap centre 
since it is a region of low surface density. This can satisfy Eq. 
m Then the inner edge of our bump function mimics the 
outer gap edge. The outer gap edge is then a potential site 
for the growth of a low-frequency one-armed spiral through 
the background torque. However, the locally isothermal re¬ 
quirement would limit this process to the outer disc, or that 
the temperature profile about the gap edge is set by the 
planet luminosity. 

Here, it is worth mentioning the transition disc 
around HD 142527, the outer parts of whic h displays 
an m = 1 asymmetry j Fukaga wa et al.l l2013h and spi¬ 
ral arms (IChristiaens et al.ll2014l 'l just outside a disc gap. 
These authors estimate Q ~1—2 in the outer disc, im¬ 
plying self-gravity is i mportant, but the disc may remain 
gravitationally-stable llChristiaens et akl l2014l l. This is a 
possible situation that our disc models represent. 


8 SUMMARY AND CONCLUSIONS 

In this paper, we have described a destabilising effect of 
adopting a fixed temperature profile to model astrophysical 
discs. By applying angular momentum conservation within 
linear theory, we showed that a forced temperature gradient 
introduces a torque on linear perturbations. We call this 
the background torque because it represents an exchange of 
angular momentum between the background disc and the 
perturbations. This offers a previously unexplored pathway 
to instability in locally isothermal discs. 

In the local approximation, we showed that this back¬ 
ground torque is negative for non-axisymmetric trailing 
waves in discs with a fixed temperature or sound-speed pro¬ 
file that decrease outwards. A negative background torque 
enforces low-frequency non-axisymmetric modes because 
they are associated with negative angular momentum. 

We demonstrated the destabilising effect of the back¬ 
ground torque by carrying out direct numerical hydrody¬ 
namic simulations of locally isothermal discs with a self- 
gravitating surface density bump. We find such systems 
are unstable to low-frequency perturbations with azimuthal 
wavenumber m = 1, which leads to the development of an 
one-armed trailing spiral that persist for at least O(IO^) or¬ 
bits. The spiral pattern speed is smaller than the local disc 


rotation and growth rates are O(10“^D) which gives a char¬ 
acteristic growth time of 0(10) orbits. 

We used three independent numerical codes — FARGO 
in 2D, ZEUS-MP and PLUTO in 3D — to show that the 
growth of one-armed spirals in our disc model is due to 
the imposed temperature gradient: growth rates increased 
linearly with the magnitude of the imposed temperature 
gradient, and one-armed spirals did not develop in strictly 
isothermal simulations. This one-armed spiral instability can 
be interpreted as an initially neutral, tightly-wound m — 1 
mode being destabilised by the background torque. The spi¬ 
ral is mostly confined between two Q-barriers in the surface 
density bump. We find the instability behaves similarly in 
2D and 3D, but in 3D the spiral disturbance becomes more 
radially global away from the midplane. 


8.1 Speculations and future work 

There are several issues that remain to be addressed in fu¬ 
ture works: 

Thermal relaxation. The locally isothermal assumption 
can be relaxed by including an energy equation with a source 
term that restores the disc temperature over a characteristic 
timescale Leiax. Preliminary FARGO simulations indicate 
a thermal relaxation timescale treiax < is needed 

for the one-armed spiral to develop. However, this value is 
likely model-dependent. For example, a longer treiax may 
be permitted with larger temperature gradients. This issue, 
together with a parameter survey, will be considered in a 
follow-up study. 

Non-linear evolution. In the deeply non-linear regime, 
the one-armed spiral may shock and deposit negative angu¬ 
lar momentum onto the background disc. The spiral ampli¬ 
tude would saturate by gaining positive angular momentum. 
However, if the temperature gradient is maintained, it may 
be possible to achieve a balance between the gain of neg¬ 
ative angular momentum through the background torque, 
and the gain of positive angular momentum through shock 
dissipation. We remark that fragmentation is unlikely be¬ 
cause the_£0;Uotation_£adiu£js_outeideHiebulk_of_the_S£ir^ 
arm JPurisen, Hartguist fc Pickettll2008l: [Rogers fc WadslevI 

[2012h . In order to study these possibilities, improved numer¬ 
ical models are needed to ensure total angular momentum 
conservation on timescales much longer than that considered 
in this paper. 

Other applications of the background torque. The back¬ 
ground torque is a generic feature in discs for which the 
temperature is set externally. It may therefore be rele¬ 
vant in other a strophysical contexts. One possibili ty is 
in Be star discs ([Rivinius. Carciofi fc MartavanI r2013[ l. for 
which one-armed oscillations may explain long-timescale 
variations in their emission lines (see e.g . [Okazakil [l997[ : 
[Papaloizou fc Savoniid [20061 : [Qgilvid [2008I . and references 
therein). These studies invoke alternative mechanisms to 
produce neutral one-armed oscillations (e.g. rotational de¬ 
formation of the star), but consider strictly isothermal discs. 
It would be interesting to explore the effect of a radial tem¬ 
perature gradient on the stability of these oscillations. 
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APPENDIX A: THE BACKGROUND TORQUE 
DENSITY IN A THREE-DIMENSIONAL DISC 
WITH A FIXED TEMPERATURE PROFILE 


We give a brief derivation of the angular momentum ex¬ 
change between linear perturbations and the background 
disc. We consider a three-dimensional disc in which the equi¬ 
librium pressure and density are related by 

p = ct{R,z)F{p), (Al) 


where F{p) is an arbitrary function of p with dimensions of 
mass per unit volume, and Cs is a prescribed function of po¬ 
sition with dimensions of velocity squared. The equilibrium 
disc satisfies 




1 

pdR 


1 ^ 

dz 


94-tot 
dR ’ 
94-tot 
dz 


(A2) 

(A3) 


Note that, in general, the equilibrium rotation H depends 
on R and 2 . 

We begin with the linearised equation of motion 
in terms of the Lagrangian displacement ^ as given by 
iLin. Papaloizou fc Kiev! (Il993l l but with an additional po¬ 
tential perturbation, 




+ 2Q,z X 


Dt 


= - VJ4-d -RR{^- VQ^) 

P p^ ^ ' 


.-V ((£ + »,) 


&p Vp Jp Vp 
p p p p 


RR (^ ■ , 

(A4) 


where D/Dt = 9* J- imQ. for perturbations with azimuthal 


dependence in the form exp {imcj)), and the S quantities de¬ 
note Eulerian pertur bations. _ _ 

As explained in iLin fc Papaloizoul ll201ll 'l , a conserva¬ 
tion law for the angular momentum of the perturbation may 
be obtained by taking the dot product between Eq. IA4I and 
(—m/2)p^*, then taking the imaginary part afterwards. The 
left hand side becomes the rate of change of angular momen¬ 
tum density. The first term on the right hand side (RHS) 
becomes 
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--Im 


m 

= -Im 


-pC ■ V 


Sp 


+ J4>d 


V- 


pe(j+ 


47rCr 


+ = In.(vf). (A5) 

where Sp — —V • (p^) and V^J4>d = AnGSp have been used. 
The terms in square brackets on the RHS of Eg. IA5l is (mi¬ 
nus) the angular momentum flux. The second term on RHS 
of Eq. IA5I together with the remaining terms on the RHS 
of Eq. lAl constitutes the background torque. That is. 


Tbg = 



<5p,* „ ... 9(RH") 


(A6) 


where Ap = Jp + ^ • Vp is the Lagrangian density perturba¬ 
tion. 

So far we have not invoked an energy equation. For 
adiabatic perturbations Tbg is zero, and we recover the 
same statement of angular m omentum conservation as in 
iLin. Papaloizou fc Kiev! lll993ll but modified by self-gravity 
in the fluxes. 

However, if we impose the equilibrium relation Eq. ED 
to hold in the perturbed state, then 

Sp= cl{R,z)F'6p, (A7) 


where F' — dF/dp. Inserting this into Eg. IA6I we obtain 


TBG = -^4lm 
2 pci 


SpC-VcI+Cr^z 


( dp del 

[di'M 


dp 9c^ Y 
dR dz J 
(A8) 


where the equilibrium equations were used. At this point 
setting = 0 gives Tbg for perturbations with no vertical 
motion, 

Tbg, 2 D = Im [SpCRdRcl) , (A9) 

Z pCg 

and is equivalent to the 2D expression, Eq. m with Sp re¬ 
placed by STj and p = c^p. 

In fact, we can simplify Eq. IA8I in the general case by 
using Sp = —pV ■ ^ ^ ■ Vp, giving 

Tbg = ^^ Im [p (V ■i)e- Vc^] . (AlO) 

Z pCs 

For a barotropic fluid p = p(p), the function cl can be taken 
as constant lEo. lAll) for which Tbg vanishes. When there is 
a forced temperature gradient, Eq. lAlOl indicates a torque 
is applied to compressible perturbations (V ■ ^ yf 0) if there 
is motion along the temperature gradient (^ ■ Vc^ 7 ^ 0). 
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APPENDIX B: RELATION BETWEEN 
HORIZONTAL LAGRANGIAN 
DISPLAGEMENTS FOR LOCAL, LOW 
FREQUENCY DISTURBANCES 


Here, we aim to relate the horizontal Lagrangian displace¬ 
ments and in the local approximation. Using the local 
solution to the Poisson equation 


|/C| 

llShulll99lll . the linearised azimuthal equation of motion be¬ 
comes 


i<jSVtprn 



27vGE\ 

~W) 


5E 


m • 


(B2) 


Next, we replace the surface density perturbation S'Em = 
—ifeE^ij, and use the expressions 


SvRm = -ifr^fl, (B3) 

Sv^m = SvRm (B4) 

a 

llPapaloizou fc Pringlel^~985^ to obtain 

- 2ian^R = Ch, (B5) 

where the dispersion relation Eq. E] was used. In the local 
approximation, \kR\ ';$> m hy assumption. Hence the RHS 
of this equation can be neglected. Then 


a 


(B 6 ) 


For low-frequency modes we have a ~ —mfl, so ~ 
‘2i^R/’m, as used in the main text. 


APPENDIX C: THE CONFINED SPIRAL AS 
AN EXTERNAL POTENTIAL 


Let us treat the one-armed spiral confined in i? G [i?i, R 2 ] 
as an external potential of the form i&ext (R) cos (0 — Qpt). 
We take 

(Cl) 

H 

where ARing is the disc mass contained within R G [Ri, R 2 ], 
R = (Ri -|- R2)I2 is the approximate radial location of 
the spiral, is the Laplace coefficient and /? = R/R. 

This form of "Fext is the m = 1 component of the gravita- 
tional potential of an extern al satellite on a circular orbit 
l|Goldreich fc Tremainel[l979l L 

We expect the external potential to exert a torque on 
the disc at the Lindblad and co-rotation resonances. At the 
outer Lindblad resonance (OLR), this torque is 


Fi, = 


TT^Sl p d^Fext 

dR 


+ 2 
L 




(C2) 


where a Keplerian disc has been assumed and subscript L de¬ 
notes evaluation at the OLR, R = Rr. (The inner Lindblad 
resonance does not exist for the pattern speeds observed in 
our simulations.) 

If we associate the external potential with an angular 


_2 

momentum magnitude of Jext = ARingR S 2 p, we can calcu¬ 
late a rate of change of angular momentum = F/^/dext- 
Then 

Qp 3 Ql\mY\r)\rJ \RcJ 

Rl dbi/2 

'W dp 

(C3) 
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Inserting h = 0.05, Ql = 10, Mrfng = 0.05M*, Rl = 7.2Ro, 
Rc = 4 . 4 R 0 and R = 1.5Ro from our fiducial FARGO sim¬ 
ulation, we get 

71 , ~ 5 X 10“'‘!2p. (C4) 


For the co-rotation torque, we use the result 


Fc = 7r"4-L 


f ^Y^ — 

VIR J dR 



(C5) 


for a Keplerian disc, where subscript c denotes evaluation at 
co-rotation radius R = Rc- For a power-law surface density 
profile E oc R”® we have 


7 c 4 7r/i 

Y ^ 3^ 



(C 6 ) 


Using the above parameter values with s = 2 and Qc = 10, 
we obtain a rate 


7c ~ 6 X 10"'‘np. 


(C7) 


The torque exerted on the disc at the OLR by an 
external potential is positive, while that at co-rotation 
depends on the gradient of potential vorticity there 
(iGoldreich fc TremaindflOLgl l. For our disc models with sur¬ 
face density E oc R~^ in the outer disc, this co-rotation 
torque is positive. This means that the one-armed spiral 
loses angular momentum by launching density waves with 
positive angular momentum at the OLR, and by applying 
a positive co-rotation torque on the disc. In principle, this 
interaction is destabilising be cause the one-armed spir al has 
negative angular momentum dLin fc Papaloizoull201ir) . 

However, the above estimates for 71 , and 7 c are much 
smaller than that due to the imposed temperature gradient 
as measured in the simulations (7 ~ O.lflp). We conclude 
that for our disc models, the Lindblad and co-rotation reso¬ 
nances have negligible effects on the growth of the one-armed 
spiral in the inner disc (but it could be important in other 
parameter regimes). 


























